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Abstract 

The energy variance extrapolation method consists in relating the ap- 
proximate energies in many-body calculations to the corresponding energy 
variances and inferring eigenvalues by extrapolating to zero variance. The 
method needs a fast evaluation of the energy variances. For many-body 
methods that expand the nuclear wave functions in terms of deformed Slater 
determinants, the best available method for the evaluation of energy vari- 
ances scales with the sixth power of the number of single-particle states. 
We propose a new method which depends on the number of single-particle 
orbits and the number of particles rather than the number of single-particle 
states. We discuss as an example the case of 4 He using the chiral N3LO 
interaction in a basis consisting up to 184 single-particle states. 
Pacs numbers: 21. 60.De, 24.10.Cn, 27.10.+h 
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1 Introduction. 



Both the Monte Carlo shell model (MCSM) (refs. [l]-[3]) and the Hybrid Multi- 
Determinant (HMD) method (refs. [4], [5]) approximate eigenstates of the nuclear 
Hamiltonian, although with a different parametrization and minimization method, 
with a linear combination of Slater determinants as a variational ansatz. Several 
years ago (refs.[6]-[8]) a robust method has been proposed whereby the error in 
the energy E approx — E exact was shown to have a linear behavior as a function of 
the energy variance for wave-functions sufficiently close to the exact ones. By ex- 
trapolating to zero variance one can then infer the exact (or almost exact) energy 
eigenvalues. This idea was put forward for the first time in the context of shell 
model calculations where the energy variance is generated by the Lanczos diago- 
nalization method. In ref. [6] a sequence of approximate shell model wave func- 
tions \ip > was obtained by truncation of number of possible excitations. Using 
the proportionality relation < ip\H\ip > —E exact oc< i/j\H 2 \iIj > — < ip\H\ip > 2 , 
valid if the approximate wave function is sufficiently close to the exact one, and 
extrapolating to zero variance one can obtain the exact value of the energy, at least 
in principle. Within the Lanczos diagonalization method the variance is obtained 
without extra effort. For variational methods, although the extrapolation method is 
still valid, the evaluation of the energy variance is computationally very expensive. 
Note however that access to the full Hilbert space is not needed. Recently the ex- 
trapolation method has been revived within the Monte Carlo shell model method 
(refs. [9], [10]) and applied to ab-initio calculations. The cost of the method is pro- 
portional to the sixth power of number of single-particle states. The very same 
cost applies to the Hybrid Multi-Determinant method. Clearly for a systematic 
applicability of the method, the computational cost must be reduced. In this work 
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we propose an efficient method to evaluate the energy variance. More precisely 
we give an efficient recipe to evaluate < ip\H 2 \(j) > where \ip > and \<p > are two 
Slater determinants. The method which will be described in detail in the next sec- 
tion is based on a factorization of the density matrix and its computational costs 
depends on the number of particles and of single-particle orbits n, /, j rather than 
on the number of single-particle states n, l,j, m. In section 3 we discuss the ap- 
plication of the method to the case of 4 ife using the N3LO interaction (ref.[l 1]) 
in an harmonic oscillator basis of 6, 7 and 8 major shells (N ho = 5, 6, 7 respec- 
tively). We restrict ourselves to I < 5 and the largest single-particle space consists 
of 184 single-particle states. Rather than work with the bare N3LO interaction we 
renormalize the interaction with two methods. In the first method we renormal- 
ize the interaction with a sharp relative momentum cutoff K max = 2.5/m~ 1 as 
done in Viowk (ref. [12]), then we renormalize once more to a specified num- 
ber of harmonic oscillator shells using the Lee-Suzuki renormalization procedure 
(refs.[13],[14]). In the second method we use directly the Lee-Suzuki procedure 
on the N3LO interaction. The method we propose to evaluate the energy variance 
has been implemented on a personal computer with modest resources. 



2 Evaluation of the energy variance. 

In both the MCSM and in the HMD methods the variational ansatz for the eigen- 
states is 

N 

IV'JV >= 9aN\U a ,N >, (1) 



a=l 
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Only for simplicity we omit the projector to good angular momentum and parity. 
The \U aN > are variational Slater determinants or, more precisely, a product of a 
neutron and a proton Slater determinant. N is the number of Slater determinants. 
If Slater determinants are added one by one the cost of generating and optimizing 
N Slater determinants scales as iV 2 . Therefore it is computationally expensive to 
increase the accuracy in the evaluation of energies and observables after a large 
number of Slater determinants. Therefore an efficient and theoretically robust 
extrapolation method is highly valuable. The coefficients g a:N are determined by 
solving the generalized eigenvalue problem 

£ < U atN \H\Uf} ,N > g/3,N — En < Ua t N\U/3 t N > <7/3,JV, (2) 

P 

E N being the approximate energy for a set of N Slater determinants. We write the 
Hamiltonian as 

H = ]^^Vi jM a\a)aia k (3) 

where a J is the creation operator in the single-particle state i. For convenience we 
include the single-particle Hamiltonian in the anti- symmetric part of the two-body 
potential v^^i = \{ v ijki — Vijik)- The matrix elements of H and of H 2 between 
two different Slater determinants | U > and | V > are given by 

<V\H\U> _ 

< y^j > = VijklPkiPlj (4) 

and 

< > = {Vijklpkipljf + 4tr( VF Vp) + VijklFkpFigVpqrsPriPsj (5) 

where p ki =< V\a\a k \U > is the generalized density matrix, F ik =< V\ai<i\\U > 
and 

Vji = Vij k ip k i (6) 



in eqs. (4)-(6), the sum over repeated indices is understood. Note that the indices 
of the density matrix refer only to identical particles, that is p np = 0. 
The use of the anti- symmetric part of the Hamiltonian v reduces the number of 
terms since the exchange terms are opposite to the direct contributions. In eq. (5), 
the trace term is taken in the single-particle indices, i.e. 



In ref.[9], apart a slight change in the notations, the last term in eq. (5) is recast as, 



and it is treated as two products of matrices of dimensions equal to N^ p , where 
N sp is the number of single-particle states. Hence the computational cost scales 



We now proceed to modify these expressions in order to reduce the compu- 
tational cost. The first two terms in eq.(5) pose no problem and we will focus 
entirely on the third term alone which we call Q. For readability we replace latin 
indices with numbers (e.g. 1 = (ni, h,ji, ra{) etc. of dimension N sp ) 

The Slater determinants | U > and | V > are identified by the single-particle 
wave functions U 1>a and V 1>a with and a = 1, 2, ..A where A is the number of 
particles (we will treat explicitly neutrons and protons later). The expressions for 
p and F are given by 



triVFVp) = VijFjkVupu 



(7) 



Q = V(i j )^ki){.FF) {k i ) ^ pq) V(p q) ,( rs )(pp) {rs)> ( ij ) 



(8) 




Pl2 = U laW a , 2 



(9a) 



and 



Fi2 



= 5u - Pl2 



(96) 



where 5 is the Kronecker 5, and 



w a , 2 = E(^^X 2 



(10) 
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From eq.(9a) and eq.(10), one can see that p 2 = p and trp = A and the matrix 
multiplications in the expression for p are of the type (AM.) x (AA) x (AN). The 
rectangular matrices U, W are stored at the beginning of the calculations. The 
various nn, pp and np contributions give 

Q = Q(nn) + Q(pp) + 4 Q(np) ^ 

the factor of 4 comes from all possible exchanges between n and p indices, and 
Qfap) contains only terms like < np\v\np >, for sake of argument, in this order. 
Inserting eqs. (9a)-(10) in the last term of eq.(5) (the Q term) and using eq.(9b), 
we obtain 

Q(tt') = v 1234 F 35 F 46 v 567 sp 7 iP82 = C 2 (tt') + C 3a (tt') + C 3b (tt') + C 4 (tt') (12) 
for (tt') = (nn), (pp), (np), with 

C 2 (W) = (W ll W S2 v 1234 )(v 34V2/ U 1 , 1 U 2/s ) (13a) 
C 3a (tt') = -(W^A¥ 82 v X2 ^U^)(W^v^ V2l U Vl U 2l& ) (136) 
C 3b (tt') = -(W^W 52 v 12U U 3a )(W a3/ v 3/4V2/ U Vy U 2/5 ) (13c) 

C A (tt') = (W r 7 iWi 2 Ul234^3a^)(W Q3 /W>4'U3'4'l'2'^l' 7 ^2'«) (13d) 

Notice that in the case (tf ) = (np), odd numbers/letters represent neutrons (e.g. 
135, .., 0:7 ..) while even numbers/letters represent protons (e.g. 2, 4, 6, ..(55, ..). 
Again the sum over repeated indices is understood. In eq.(13a)-(13d) we grouped 
the various terms so that matrix multiplications can be efficiently performed. No- 
tice that 

£ 7 5,34 = W yl Ws 2 Vl 23 4 (14a) 

and 

-^34,7.5 = V 3 4y 2 >Ui>jU 2 >5 (146) 
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enter in all contributions and can be very efficiently evaluated since v is very 
sparse. The various C 2 , C 3 , C 4 carry a label that identifies their 2-body, 3-body 
and 4-body character. Each term contained between the brackets can be evaluated 
from the previous equations and using eqs.(14a),(14b). Eqs. (11)-(14) improve 
the computation of the Q term, since some of the indices are particle indices, 
rather than single-particle ones. Essentially we make use of the fact that p is a low 
rank matrix. In order to see how the computational cost scales with the particle 
number and with the number of single-particle states, consider for example the 
case (tt 1 ) = (nn) and let N n be the number of neutrons and N v the number of 
non-zero matrix elements of v. It is easy to see that the computational cost for the 
evaluation once for all of the matrices L and R in eq.(14) scales as N V N% and the 
evaluation of C 2 , C 3 and C 4 scale as N^N^ N%N% p and N*N sp respectively. Let 
us remark that in ab-initio calculations N sp » N n . 

Despite this improvement, we prefer to use the angular momentum coupled 
matrix elements of v. The reason is two-fold. First of all for large single-particle 
spaces the number of uncoupled matrix elements of v is very large. In some of 
the calculations described in the next section we have used a personal computer 
of lGyb of RAM, Second, instead of dealing as in eqs.(13a)-(13d) with single- 
particle indices, we prefer to deal with orbits (i.e. nlj). 

Let us first change slightly the notations: let us label orbits by numbers, i.e. 1 = 
(ni,hiji) etc., while single-particle indices are now represented as (l,mi) etc., 
and, as before, greek letters count particles. It is then straightforward to arrive 
at the following results. Define the angular momentum coupled quantities by 
Clebsh-Gordan coefficients 

R J A = v J ul2 (UU) J Z (16) 



7 



L J t& = (WW)ft4 (17) 

with 

(WW) J a f 12 = E <jim 1 j 2 m 2 \JM>W a , lmi W p , 2ma (18) 

(UU){^= J2 < h™ij 2 m 2 \JM > U lmi , a U 2m2 ^ (19) 

mitii2 

We stress again that in the case of the (np) contribution, odd numbers/letters refer 
to neutrons and even numbers/letters refer to protons. Then the 2-body contribu- 
tion to Q is given by 

C2(U') = ^L J a f 12 R{Z (20) 

JM 

where we have shown explicitly the sum over JM only, and the sum over the 
remaining indices is implicit. All matrix multiplications have now small dimen- 
sions. For the contribution of 3 -body type, define also 

W^ M (j 3 m 3 ) = E < hmshme | JM > W Pfim<j (21a) 

m 6 

U^ M {j 3 m 3 ) =J2< J3m 3 j 4 m A \JM > U 4mi ,p (216) 



TO4 



and 



^, 75 (3,m 3 ) = £ Wff(j3,m3)R{™ s (22a) 



JM,6 



£ T ^(3,m 3 ) = £ L J ^ u Ulf(j 3 m 3 ) (226) 



JMA 



Then 



Csaitt') = -£ C l5 ^3m 3 )n^m 3 ) (23) 

3m 3 

In the case of C 3b define 

Ui™{jim 4 ) = J2< hm 3 j4m 4 \JM > U 3ms>a (24a) 
(j 4 m 4 ) = E < j 5 m 5 j 4 m 4 1 JM > W a , 5ms (246) 
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W15 



and 

£ 75iQ (4, m 4 ) = £ L™ 4 Ui M (j 4 m 4 ) (25a) 
^(4, m 4 ) = £ O4, m 4 )i? 5 J 4 M 75 (256) 

hJM 

Then 

C3 6 («') = - E ^7<5,«^,7<5( 4m 4) (26) 
4,m4 

Note the difference in the implicitly summed indices in this expression compared 
with the ones in eq.(23). The contribution of 4-body type has the simple expres- 
sion 

C 4 (ft') = M a( s r/S M^ (27) 

with 

M a ^s = T,(WW) J a f 12 R{ 2 M , S (28) 

JM 

In the np case the contributions C 3a and C 3b are not equal in general. They are 
equal in the case of (tt') = (nn), (pp). Also it is worth to note that the four-fold 
sum in eq.(28) is over the particle indices only. As before, let us analyze how the 
computational cost scales the number of orbits and the number of particles. Let 
us consider the case of the nn contribution. The evaluation of eqs.(18) and (19) is 
straightforward because of the angular momentum selection rules. The evaluation 
of eqs.(16) and (17) scales as N JM N V (JM)N% where N JM is the number of the 
possible JM values, N V (JM) is the number of non-zero matrix elements wf 234 
(the total number of angular momentum coupled matrix elementsU^ can be a 
factor of 30 smaller than the number of uncoupled matrix elements of v). Once 
eqs.(16) and (17) have been evaluated for all JM values, the computational cost of 
C 2 scales as NjmN^N^ where N orb is the number of single-particle orbits. The 
evaluation of C 3 (cf. eq.(23)) scales as N sp N^ in addition to the cost of evaluation 
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of eqs.(22a) and (22b) which scales as N JM N sp N orb N^ (note that the 3-body term 
retains one power of N sp . The computational cost of the evaluation of the matrix 
M. scales as N JM Nl rb N^ and the evaluation of the 4-body contribution scales as 
N*. Since the number of uncoupled matrix elements of the Hamiltonian can easily 
be of the order of 10 7 we use only the angular momentum coupled representation. 



3 A case study with A He. 



Let us apply the method derived in the previous section to 4 He. We use the 
N3LO interaction (ref. [11]) as the nucleon-nucleon potential. This potential 
has a smooth energy cutoff at A = 500MeV corresponding to a relative momen- 
tum K max = 2.534/m~ 1 The potential is however non-zero at higher relative 
momenta due to the smooth cutoff. In ref. [15], a rather large number of harmonic 
oscillator major shells was used in order to reach independence from the h.o. fre- 
quency and from the single-particle space, using the 'bare' potential. In order to 
be able to work with smaller spaces we renormalize the interaction with two meth- 
ods. In the first method we first renormalize the interaction much in the same way 
it is done to obtain Vi owk interactions (cf. ref. [12]) to a sharp relative momentum 
K max = 2.5/m -1 (essentially to get rid of the high momentum tail) and then we 
renormalize once more to a specified number of major shells with the Lee-Suzuki 
method, in order to reach some independence from the h.o. frequency with rea- 
sonably small single-particle spaces. In the second method we apply directly the 
Lee-Suzuki method to the N3LO interaction. 
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The rationale for the first choice is the following. For an interaction which has 
a sharp cutoff K max and for a system of radius R the acceptable values of Ml are 
given by the inequality (ref.[15]) 

tfK 2 m J{mN ho ) <hu< N ho h 2 /(mR 2 ) 

and in order to have a reasonable large interval in kui we need large values of N ho . 
This is correct if we use the bare interaction. It is hoped that, by renormalizing the 
interaction one more time, we can reach the energies given by the bare interaction 
having a specified K max , with smaller values of N ho . The only problem with 
this argument is that the first renormalization to a sharp cutoff could give rise to 
induced many-body interactions. This is precisely the reason why we took a large 
value of K max , that is, to minimize the induced many-body interactions. In ref. 
[17], it was shown that in order to minimize the effect of induced many-body 
interactions, K max should be as large as possible (K max ~ 5/m" 1 ), this result 
has been obtained using the CD-Bonn interaction which contains larger energy 
scales than N3LO. In any case, if we use Vi ow k interactions with a large cutoff, 
another renormalization step is nearly unavoidable, unless we are willing to work 
with very large single-particle spaces. 

We have used the HMD method of type-a (described in ref. [16]) to generate 
a sequence of wave-functions \ip n > consisting of a n Slater determinants (up to a 
maximum of 100). We have used a projector to good parity and good z-component 
of the angular momentum. We considered several values of the harmonic oscil- 
lator frequency namely htt = (21, 24, 27, 30, 33, 36)MeV and N ho = 5, 6. For 
Nho = 7 we considered only Ml = (27, 30)MeV. Only states with I < 6 have 
been included in the single-particle basis. The largest number of single-particle 
states is 184 and the corresponding number of orbits is 32 for N ho = 7. 
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Figure 1 : Energy-variance plot for N ho = 7. 



12 



-22 



-28 











'V lowk+LS N ho=5' 












'V lowk+LS N ho=6' 


— X— 










'V lowk+LS N ho=7' 


— e— 










IS N ho=5' 


b 


- 








'FY' 






























X— 






. c ' 












i i 





20 25 30 35 40 

*!2 (MeV) 

Figure 2: Energies as a function of Ml 

Among the several variants of the energy vs. variance extrapolation methods 
we considered the first one, used in ref.[6], where E n =< ip n \H\ip n > is expanded 
as a function of the dimensionless variance defined as 

AE = Q 2 = < jj n \H 2 \lp n > - < 1p n \H\jj n > 2 

< 1p n \H\ljj n > 2 < ^ n \H\lp n > 2 

where \ip n > is the wave-function obtained with the first n Slater determinants, 
and fitted with a linear function E n = a + bAE n . In the limit AE n = 0, a is the 
extrapolated ground-state energy. 

For the purpose of the extrapolation we discard low values of n as done in ref. 
[9]. As a measure of the soundness of the linear fit we consider the quantity 

S no = iY,{E n -a- bAE n ) 2 (29) 

y n>no 
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where the first n many-body wave functions have been discarded, since they are 
a poor approximation to the exact eigenstate. The value of n is selected so that 
£ n0 is of the order of lOKeV or smaller. In fig. 1 we show a few plots of the 
energy as a function of the variance a 2 . Typically variances for this nucleus and 
this interaction are much larger than the ones shown in ref. [9]. This is not very 
surprising since in our case the spectrum extends over very large values of the 
energy, while in the case study of ref. [9], where only one and two major shells 
were considered, the spectrum is much more compressed. In fig. 2 we show our 
results of the calculation and the the exact Faddeev-Yakubovski result. In fig. 2 we 
also show the results obtained using directly the Lee-Suzuki method for N ho = 5. 
Interestingly enough, the double renormalization method performs better than the 
single L-S renormalization which shows a too marked variation as a function of 
the harmonic oscillator frequency. Note also the results for N ho = 7 are very close 
to the ones for N ho = 6. 

In conclusion, in this work we have implemented an efficient method to eval- 
uate energy variances for variational many-body calculations which are needed in 
order to apply the energy variance extrapolation method. This method avoids the 
need to determine a very large number of Slater determinants in order to improve 
the accuracy of observables. Since the computational cost of the evaluation of N 
Slater determinants scales as N 2 , the method is a useful tool to decrease the cost 
of variational many-body calculations. 
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